Exact dynamics of interacting qubits in a thermal environment: 
Results beyond the weak coupling limit 
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We demonstrate an exact mapping of a class of models of two interacting qubits in thermal reservoirs to two 
separate spin-bath problems. Based on this mapping, exact numerical simulations of the qubits dynamics can 
be performed, beyond the weak system-bath coupling limit. Given the time evolution of the system, we study, 
in a numerically exact way, the dynamics of entanglement between pair of qubits immersed in boson thermal 
baths, showing a rich phenomenology, including an intermediate oscillatory behavior, the entanglement sudden 
birth, sudden death, and revival. We find that stationary entanglement develops between the qubits due to their 
coupling to a thermal environment, unlike the isolated qubits case in which the entanglement oscillates. We also 
show that the occurrence of entanglement sudden death in this model depends on the portion of the zero and 
double excitation states in the subsystem initial state. In the long-time limit, analytic expressions are presented 
at weak system-bath coupling, for a range of relevant qubit parameters. 
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I. INTRODUCTION 

Understanding the dynamics of a dissipative quantum sys- 
tem is a prominent challenge in physics, as a quantum sys- 
tem is never perfectly isolated from a larger environment. A 
minimal, yet highly rich model for exploring quantum dissi- 
pation effects, is the spin-boson model, including an impurity 
two-level system (referred to as a spin) coupled to a thermal 
reservoir. This model displays a rich phase diagram in the 
equilibrium regime U Q] • The non-equilibrium version of 
this model, referring to the case where the spin is coupled 
to two thermal reservoirs, has been suggested as a toy model 
for exploring quantum transport phenomenology through an 
anharmonic nanojunction In this case, the generic sit- 

uation is one of a non-equilibrium steady-state, regardless of 
the initial preparation. 

Interacting two-level systems are the basic element in quan- 
tum computation, thus it is paramount to extent the minimal 
spin-boson scenario and describe more complex modular sys- 
tems, e.g., two interacting qubits immersed in a thermal envi- 
ronment J5t]. For a schematic representation, see Fig. Q] The 
qubits may share their thermal environment, or may separately 
couple to independent baths, maintained at a nonzero temper- 
ature. The latter situation corresponds to the case where the 
qubits are not necessarily placed close to each other. In an- 
other relevant setup, one qubit couples indirectly to a thermal 
reservoir, through its interaction with the other qubit. This 
situation effectively corresponds to a subsystem anharmoni- 
cally coupled to a harmonic bath, allowing to introduce non- 
trivial nonlinear effects [6]. Physical realizations include, for 
example, ultracold atoms in optical lattices B7D, trapped ions 
Jit], resonator-coupled superconducting qubit arrays l9l [Toll , 
and electron spins in quantum dots and doped semiconduc- 
tors 11 ill . In such systems one should consider (at least) 
four energy scales: the internal qubit energetics, controlling 
its isolated (Rabi oscillation) dynamics, qubit-bath interaction 



strength and the environment temperature, leading to decoher- 
ence and relaxation processes, and qubit-qubit coupling en- 
ergy, admitting state transfer between qubits and a nontrivial 
gate-functionality. 

The dissipative multi-qubit system has recently served as 
a simple model for resolving issues related to coherence dy- 
namics in the time evolution of biological molecules, e.g., the 
Fenna- Matthews-Olson (FMO) complexes, resulting in the 
identification of the relevant decoherence, relaxation and dis- 
entanglement timescales Ill2l4l4ll . It should be noted that these 
works have considered single-excitation states only, ignoring 
the contribution of the zero and the doubly excited states in 
the two-qubits dynamics. 

In this paper, we analytically demonstrate that a class of 
interacting two-qubit systems immersed in separate thermal 
reservoirs or within a common bath, can be mapped onto two 
uncoupled spin-bath problems, allowing for an exact numer- 
ical solution of the qubits dynamics. For a bosonic environ- 
ment and a particular system-bath interaction form, we per- 
form those simulations using an exact numerical technique, 
the quasi adiabatic path-integral (QUAPI) approach IU5I1 . pro- 
viding the population and coherence dynamics of the system. 
With this at hand, we can follow the exact dynamics of entan- 
glement between the qubits, as quantified by Wootters' con- 
currence 11611 . For a range of system and bath parameters, 
analytic results are presented, describing the system behavior 
in the long time limit. The model investigated here is more 
general than what has been typically considered before, going 
beyond the simple exchange interaction model lfl3l[l4ll . 

Entanglement is associated with nonclassical correlations 
between two or more quantum systems [ 17]. Since it is a ba- 
sic resource in quantum computation and information tech- 
nology, it is important to understand the extent to which 
environmental-induced decoherence processes degrade and 
destroy it lfl8ll . or alternatively, generate lfl9tl and maintain 
it ll20ll . It has been recently shown that two qubits in separate 
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FIG. 1 : Scheme of the model system including two interacting qubits 
(a) immersed in a common bath, (b) coupled to separate baths, L and 
R, and (c) with qubit '2' coupled to a thermal bath only through its 
interaction with qubit '1'. This model can represent the nonlinear 
coupling of qubit '2' to a structured bath. Simulations were per- 
formed here assuming scenario (c). 



reservoirs may disentangle at finite times, as opposed to the 
behavior of coherences. This process is referred to as "entan- 
glement sudden death" lfl8il2lll . More recent theoretical and 
experimental studies have looked at related effects, e .g., the 
collapse and subsequent revival of the entanglement II22(|. or 
its delayed-sudden birth, induced by a dissipative bath B23I1 . 
Steady-state entanglement generation by dissipation has been 
recently observed in atomic ensembles [24]. These studies 
have assumed non-interacting qubits, and the system dynam- 
ics has been typically followed within quantum master equa- 
tion approaches (e.g., the Redfield equation or Lindblad for- 
malism ||25ll ), by invoking the weak system-bath coupling ap- 
proximation. The markovian limit has been further assumed 
in many cases, see e.g., II191 12611 . 

Our work here departs from these studies in two substan- 
tial aspects. First, we consider a more complex model for the 
subsystem, introducing qubit-qubit interaction, with the mo- 
tivation to examine a setup relevant for quantum computing 
technology. Using an exact mapping, we show how the dy- 
namics can be followed within a simpler construction: While 
we take into account the zero excitation and double exci- 
tation states, under certain initial conditions their dynamics 
can be separated from the evolution of the single-excitation 
states. However, their contribution to the pair entanglement 
is paramount. Second, we refrain from making approxima- 
tions: We study the qubits dynamics using a numerically ex- 
act method, assuming a class of initial states. The results are 
valid beyond the weak system-bath coupling scenario, accom- 
modating non-markovian effects. 

Our calculations display rich dynamics. Particularly, we 
observe the development of a stationary concurrence due to 
the coupling of the qubits to thermal baths. This result stands 
in a direct contrast to the oscillatory behavior observed in the 
fully coherent regime. It demonstrates that while entangle- 
ment inherently relies on the existence of quantum correla- 
tions in the system, it nevertheless requires non-vanishing de- 
coherence and relaxation effects in order to be stabilized and 
become useful for quantum technologies. Other phenomena 
detected and explained here are entanglement delayed sudden 



birth, sudden death, and revival. 

The paper is organized as follows. In Sec. II we describe 
the model of interest, and explain its mapping onto two sep- 
arate spin-bath problems. In Sec. Ill we explain how we fol- 
low the system dynamics, and include relevant expressions for 
calculating the qubits concurrence. Numerical results within 
QUAPI are included in Sec. IV. The long-time limit is dis- 
cussed in Sec. V. Sec. VI concludes. 



II. MODEL 

The general model to be considered here includes two inter- 
acting qubits, i = 1,2, immersed within separate reservoirs, 
L and R, respectively. The formalism can be reduced to de- 
scribe a single-bath scenario. The total Hamiltonian includes 
three terms, 

H = Hs + H b + V S b. (1) 

Hs and Hb = Hl + Hr stand for the system and reservoirs 
Hamiltonians, respectively. The former includes the isolated 
qubits with the internal energy bias and a qubit-qubit inter- 
action term V ss , 

H s = eiof + £2^2 + V as , 

Vss = ^[(l+7)(7M + (l-7VM + 2Ma|]. (2) 

of (p — x, y, z) are the Pauli matrices for the ith spin, J is 
an energy parameter characterizing the exchange interaction, 
7 and (5 set the interaction anisotropy. Our mapping holds for 
a dephasing-type system-bath interaction model, 

V SB = afB L + a%B R . (3) 

Here, B v is a v = L,R bath operator, with Bl coupled to 
spin '1' and Br coupled to spin '2'. In our simulations below 
we adopt bosonic reservoirs: each thermal baths includes a 
collection of independent harmonic oscillators, 

H u = ^2uj k a\a k . (4) 

The operators a\ and a k are bosonic creation and annihilation 
operators, respectively, Uk is the mode frequency. We also 
assume that the interaction operators constitute the reservoirs 
displacements from equilibrium, 

B - = Yl Xk (4 + «*) ■ < 5 > 

kei' 

Here, are system-bath coupling constants. The mapping 
described next, from the Hamiltonian (Q]i into two spin-bath 
problems, neither rely on a particular bath statistics nor on 
the details of the operators B v . For example, it is valid for a 
model of two-qubits in fermionic or spin environments. The 
Hamiltonian introduced so far takes into account two inde- 
pendent reservoirs. We could explore a similar setup with one 
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qubit coupled to a thermal reservoir indirectly, through its in- 
teraction with the second qubit, mimicking nonlinear effects 
jfj, see Fig. [T] Another relevant setup includes two qubits 
immersed in a common thermal reservoir. 

The Hilbert space of the qubits is spanned by four vec- 
tors, 1 00) , 1 01} , 1 10} , |11), forming the excitation basis of the 
Hamiltonian, where the left (right) digit indicates the state of 
qubit 1 (2). We now show that the Hamiltonian can be mapped 
onto two spin-bath type models. We begin by defining four 
composite system operators, 

Pz = ~W-*f), Q. = \(o*i+°Z)> 

Pec = \ (of at + a\al) , Q x = \ (afa% - afa y 2 ) . (6) 
In the excitation basis, these operators take the explicit form 

P x = |01) (01| + |10) (10|, P z = [10) (01| - [01) (10| 
Qx = |00} (11| + |11) (00|, Q z = |00} (00| - |11) (11|. (7) 

Additionally, we construct two identity-type operators, 

Ip = |10)(01| + |01)(10| 

I Q = |00)(00| + |11)(11|. (8) 

With these at hand, the Hamiltonian (HJ can be written as 



H = Hr 



(9) 



where 



H Q = eQ z + JjQ x + Q Z B + J5I Q + H B I Q 

Hp = IP Z + JP X + P Z B - J5I P + H B I P . (10) 

Here, e = ei + €2, e = t\ — € 2 , B = Bl + Br and 
B = Bp — Br. One can easily show that the following com- 
mutators vanish [27] (m, n — x,z) 



[Qm,Pn] = [Pm,<Tt<T$\ = [Qm,^!] = 0. 



(11) 



The Hilbert space of two qubits can thus be factored into two 
direct-sum subspaces: The first, P, is spanned by 1 01} and 
1 10). The second, Q, is spanned by |00) and|ll). One can 
further prove that P x , [P X ,P Z ] and P z generate an SU P (2) 
group, and Q x , [Q X ,Q Z ] and Q z generate another SU®(2) 
group. The two groups have a direct-sum structure, SU P (2)Q) 
SU®(2). We now note that P x and P z in subspace P play the 
role of the Pauli matrices a x and a z (in the space spanned by 
|0) and |1)). The same principle holds for Q x and Q z in the 
Q subspace. Overall, in mapping Eq. (Q~|) into Eq. dTob we 
replaced a model of two interacting spins coupled each to its 
own thermal reservoir, by a model of two separate spin-boson- 
type systems, where each spin in the new model is coupled 
to both reservoirs. The latter model is significantly simpler 
than the former, and we can explore its dynamics using exact 
simulation tools for a class of certain initial conditions. 

The mapping described here holds for general reservoirs 
and a bilinear system-bath interaction form, with an arbi- 
trary bath operator coupled to the subsystem. The results 



also hold when we apply a dressing transformation B28I1 W, 
H' = W^HW on the two qubits Hamiltonian, or the reser- 
voirs. For example, we may introduce a spin-orbital coupling 
into V ss via W — exp(i^af), while keeping other terms in 
the total Hamiltonian unchanged 



V'. = cos 6V S! 



J, 



sin^[(l + 7 )<7M - (1 - 7K^]-(12) 



Another relevant case that can be simulated exactly relies on 
the absence of external fields, t\ = €2 = 0, taking W = 
exp(ij(af + Orf)). This results in 

H' s = ^[(1 + l)o\ol + (1 - 7 )aM + 26afa% ] (13) 

with VL B = XicrfBi + ^a^Bp. After this transformation, 
the model has turned into the anisotropic XYZ-type model 
with flip-flop (a x ) coupling between the system and reser- 
voirs. In this form, the model describes energy exchange 
between the qubits and the baths, unlike the original Hamil- 
tonian [Eq. ([3])] which delineates dephasing effects. When 
1 — 7 = 2(5, it reduces to the standard XY model, H' s = 
^{alal + ala%) + (1 + 7)^]. 



III. DYNAMICS AND QUANTUM ENTANGLEMENT 

We explain here how we time-evolve the reduced density 
matrix, to obtain the qubits dynamics. As an initial con- 
dition we consider a system-bath product state, where the 
reservoirs are maintained in a canonical-thermal state, p v = 
e~P" Hu /Z v , Z v = TiB[e~ l3 " H ''] is the partition function. In 
order to separate the dynamics into the Q and P branches, 
we must adopt a direct sum Q-P initial state for the qubits. 
Overall, the total density matrix at time t — is written as 



p(o)=r p (0) 



pq(0) 



Pb- 



(14) 



For a particular example, see Eq. ( |23l . Under this construc- 
tion, the reduced density matrix follows 

Ps (t) = T lB [U(t)p(0)U\t)] 

' U P {t)pp{0)p B U P (t) \ 

U Q (t) PQ (0)p B U Q (t) J 

(15) 



The trace is performed over the L and R degrees of freedom. 
The time evolution operators, U(t) = Uq{£) © Up(t), are 
defined as 



u Q (t) 



-UHq 



U P {t) 



itHp 



(16) 



with Hp and Hq given in Eq. ( TTOb . Eq. (TT~5T > establishes 
an important result: The dynamics of the two-qubit system 
proceeds in two independent branches, each equivalent to a 
spin-bath model. In the case of bosonic baths, the dynamics 
in each branch is followed next using the QUAPI technique 
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01511 . a numerically exact simulation tool that can be easily ex- 
tended to include more than one thermal reservoir. The output 
of this calculation is the reduced density matrix of the qubits. 
We use this information and investigate the time evolution of 
entanglement between the qubits. As a side comment, we note 
that if the qubits were to couple to a common bath, B = 0, 
only the Q subspace would have become susceptible to deco- 
herring effects, while the P subspace would be an invariant 
subspace, or a "decoherence free" subspace. 

Based on Eqs. ([T4li-(fT6]>, we conclude that the re- 
duced density matrix p$ is an X-type matrix at all times, 
once we organize it in the standard order of basis vectors 
|00> , |01> , 1 10) and |11), 
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This is an important result: The dynamics of a class of dissi- 
pative interacting qubits [Eqs. (Q]l-([3]l] can be reached via the 
solution of two spin-bath problems, and the reduced density 
matrix satisfies an X-form at all times. Different quantities 
are of interest, ej* the timescale for maintaining coherences 
in the system 11411 . Here, we focus on quantum correlations 
in the system 11711 . computed next in a numerically exact way 
beyond weak coupling. In particular, we quantify the degree 
of entanglement between the qubits using Wootters' concur- 
rence 1160 . For mixed states it is calculated by considering the 



eigenvalues of the matrix r(t) = pg(t)af < 
given here by 



lp* s (t)al®<7, 



Al,2 — 
A3, 4 = 



J (as)oi,oi(as)io,io ± |(ps)oi,io| 
J (ps)oo,oo(ps)ii,n ± I (as)oo,ii I 



(18) 



In terms of these eigenvalues, the concurrence is defined as 



C(t) — max(2 max(y Ai, y A2, y A3, y A4), 

y r \x-y%-y%-y%,Q). (19) 

It varies from C = for a disentangled state to C = 1 for a 
maximally entangled state. In the present case it reduces to 



C(t) =max(0,2Fi,2F 2 ) 



(20) 



with 



Fi = |(ps)oi,io| - [(ps)oo,oo(ps)ii,ii] 1/2 > 

F2 = |(ps)oo,n| - [(as)oi,oi(as) 10,10] 1/2 • (21) 

The dynamics of concurrence for an X-state density matrix 
has been examined in different works. For example, in Ref. 
rf26tl it is demonstrated that the effect of entanglement sud- 
den death should always take place in a noninteracting qubit 
system once coupled to a finite temperature reservoir. As we 



mention in the introductory section, we depart from this study 
and similar works in two aspects: (i) We build the reduced 
density matrix using an exact numerical treatment, and (ii) 
we consider a more general model, including quit-qubit inter- 
action effects, with the motivation to consider a setup more 
relevant for quantum computation technologies. 



IV. NUMERICAL RESULTS 

We simulate the spin-boson dynamics in the P and Q 
branches (separately) using QUAPI 01511 . to obtain the pop- 
ulation and coherences in each branch. With this at hand, 
we generate the 4x4 reduced density matrix ps(t), Eq. 
([T7l >. The qubits degree of entanglement is calculated us- 
ing Eq. (l20b . Our general description assumes two thermal 
reservoirs: Hl, coupled to spin 1 and Hr, coupled to spin 
2. These reservoirs are characterized by the spectral func- 



(17) tion J v (u) =7rJ2 kei/ X 2 k S(e k 



oj). Specifically, we simu- 
late Ohmic baths, J v {lS) = ^we""/" 1 ; uj c is the cutoff 
frequency. The dimensionless prefactor K u is referred to as 
the Kondo parameter, describing the strength of the system- 
bath interaction energy. In practice, our simulations were per- 
formed without the R reservoir, by taking Kr — 0. The rea- 
son for this choice is that in the spin-boson Hamiltonian ( [Tol l 
the inclusion of identical reservoirs which are interacting in 
the same manner with the spins (same functional form for Bl 
and Br, up to a sign), simply amounts to an additive opera- 
tion, reflecting a linear scaling of the Kondo parameter when 
more than one reservoir is incorporated. In what follows, we 
thus use the short notation K = Kl, Kr \ = and T = Tr. 

The energy parameters in the system are the qubit-qubit 
coupling, taken as J = 1, and the anisotropy parameters 
S = 0.1, 7 = 0.5. The qubits are assumed identical with 
(\ = € 2 ~ 0.1 — 0.5. For the reservoir we take as a cut- 
off frequency lo c — 7.5, and use temperatures at the range 
T = 0.1 — 1. The Kondo parameter extends from the weak 
coupling limit (K = 0.05) to the strong-intermediate regime, 
K = 0.8, where convergence of QUAPI can be achieved. The 
following initial condition is utilized for the qubits subsystem 



p s (0) = a PQ (0) © (1 - a)Ap(0), 



(22) 



with < a < 1 and pq,p(0) as (maximally entangled) Bell 
states, 



pq(o) 



1 1 
1 1 



Pp(0) 



1 1 
1 1 



(23) 



The concurrence (|20| > can be simplified if the following con- 
ditions are simultaneously satisfied: (i) a < 1/2, and (ii) t\ = 
£2- The latter condition, combined with the initial state ascrib- 
ing identical weight to diagonal elements in the P subspace, 
implies that the populations in the P (single excitation) sub- 
space are identical at all times, (ps)oi,oi — (ps)io,io = 
Since ^/((Os)oo,oo(/°s)imi < § at a ll times and the density 



matrix positivity condition demands that \pi 



< 



Pi,iPj,. 



the off diagonal terms are bounded by |(ps)n,oo| < §• This 
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implies that F2 cannot be larger than zero at any instant for 
a < 1/2. The concurrence can then be simplified to 



C o< i(t)=max(0,2Ji) 



(24) 



This expression indicates that C is nonzero if the magnitude 
of the coherence in the P subspace is large, in comparison 
to the product of populations in the Q subspace. Thus, to 
understand the behavior of entanglement between the qubits 
one needs to follow the population and coherence dynamics 
in both subspaces. At the special point a = the concurrence 
reduces to the simple form 



C=o(i) 



(|(a 



s)oi,w 



,0 



(25) 



which only depends on coherences behavior, a continuous 
function. As a result, concurrence sudden death is eliminated, 
indicating that this effect is directly linked to the inclusion of 
zero and double excitation components in the dynamics. 
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FIG. 2: Population dynamics (a) in the single excitation subspace P 
and (b) in the zero and double excitation subspace Q. We use the 
Bell states M2\ as the initial density matrix with a — 1/2. 7^=0.05, 
0.1, 0.2, 0.3, 0.4, 0.5, bottom to top; the data at K=0.05 is the most 
oscillatory. Other parameters are ei = £2 = 0.2, 5 = 0.1, 7 = 0.5, 
T = 0.2 and J = 1. QUAPI was used with a time step 5t = 0.25 
and a memory time r c = 9St. Convergence was verified by studying 
the behavior at different time step St and memory size r c . 

The qubits population behavior in time is displayed in Fig. 
12 using a = |. The qubits have the same energy gap, thus in 
the P subspace the two states are degenerate and their popula- 
tion is identical at all times, independently of K. In contrast, 
in the Q subspace the energy difference between the states is 
significant, larger than the temperature, T/2e < 1; the tun- 
neling element is given by 7 J = |, with 7 as the anisotropy 
in the qubit-qubit coupling. In such a situation we expect the 
steady-state population of the spin-up state to be significantly 
smaller than the ground state population, as indeed we observe 
in Fig. Efb). An interesting observation is the phenomenon of 
population inversion between the zero and the double excita- 
tion states before steady-state sets in. This behavior occurs 
roughly up to a timescale that is inversely proportional to the 
Kondo parameter K, independent of the temperature. For the 
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FIG. 3: Real and imaginary parts of the coherences in the single exci- 
tation subspace P, (a) and (c), and in the zero and double excitation 
subspace Q, (b) and (d). The different lines were calculated with 
Ji"=0.05, 0.1, 0.2, 0.3, 0.4, 0.5. Parameters are the same as in Fig. |2] 
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FIG. 4: Concurrence between the two qubits as a function of time, 
manifesting a steady-state bath-induced entanglement generation. 
The different lines were calculated with A'=0.05, 0.1, 0.2, 0.3, 0.4, 
0.5. Parameters are the same as in Fig. [2] Right panel: The concur- 
rence dynamics in the absence of a thermal environment for the same 
set of qubits parameters. 



same set of parameters Fig. [3]presents the coherence dynam- 
ics in the two subspaces. Generally, coherences are diminish- 
ing with the increase of K. Given the population and coher- 
ence dynamics, we display in Fig. [4] the concurrence, calcu- 
lated using Eq. (l24"i i. manifesting a rich dynamics. The follow- 
ing characteristic's are of particular interest: (i) The birth-time 
of the concurrence, (ii) its oscillations, (iii) the occurrence of 
sudden death and revival, and (iv) the steady-state value. We 
now explain those properties. 

Time-zero concurrence. The particular initial condition 
used here, a = i, results in C(t = 0) = 0. This is because 
while we are using maximally entangled states within each 
subspace as an initial condition, the entanglement between 
the two qubits themselves is zero initially, since all relevant 
reduced density matrix elements, necessary for evaluating Eq. 
( 1241 are identical JH. 



Delayed sudden birth. When (ps) 



00,00 



(P. 



a sit- 



uation taking place at, and close to, the initial time, the con- 
currence should be zero, given the positivity condition that 
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limits the value of off-diagonal elements. For small K, the 
time it takes the system to depart from its initial-equal popu- 
lation state is prolonged compared to a large-X case, thus, the 
concurrence birth-time is delayed with respect to the large K 
behavior. Interestingly, the delay in the birth time does not ex- 
tend linearly with K. Rather, the delay is significant for both 
K = 0.05 (weak system-bath coupling) and for K = 0.5 (in- 
termediate coupling), while it is shorter for in-between values, 
K ~ 0.3; The reason is that the delay time is a nontrivial func- 
tion of both the time it takes the coherences to establish, and 
the time it takes the population to significantly depart from the 
initial (equally populated) setup. 

Oscillations. The oscillatory nature of C in time, best man- 
ifested for K < 0.1 reflects the Rabi-type oscillations of the 
diagonal elements (ps)oo,oo an d (ps)u,n- When these ele- 
ments are similar in value, the concurrence drops, and even 
dies during a certain time interval, depending on the magni- 
tude of the coherence at that time. 

Steady-state value. If the two qubits are isolated from ther- 
mal effects (K=Q, right panel of Fig |4)l, the concurrence os- 
cillations reflect the nature of the population and coherences 
dynamics, depicting Rabi oscillations. The qubits behavior 
under a dissipative thermal bath is notably distinct: Since both 
population and coherences approach a constant at long time, 
the concurrence reaches a steady-state value as well. It pre- 
dominantly reflects the magnitude of the coherence (ps)io,oi 
in the long time limit since the population weakly depends on 
K at long time, see Fig. |2jb). Interestingly, for the present 
a = | case the steady-state value of the concurrence is almost 
identical at weak system-bath coupling, K = 0.05 — 0.3. It 
significantly degrades around K = 0.4 — 0.5. Beyond that, it 
is identically zero. 

Sudden death and revival. Based on Figs. [20 we can 
draw general conclusions regarding the process of entangle- 
ment sudden death. The effect is directly linked to the exis- 
tence of population in the zero and double excitation states. 
If the dynamics within the Q space is eliminated all together 
(a = 0), the concurrence is only controlled by the absolute 
value of the coherence |(/Os)oi.io|> Eq. d2~5l l. This quantity 
does not manifest an oscillatory behavior: Under the present 
initial condition it starts at a large value, touches zero at a par- 
ticular time, then grows again to a certain extent. (Under a dif- 
ferent initial condition, e.g., |(ps)oi,io(0)| = 0, the entangle- 
ment will systematically grow, up to the steady-state value). 
In contrast, when the double excitation state is initially popu- 
lated, oscillations between the states in the Q subspace largely 
occur, if system-bath coupling is weak. Then, the competition 
between the two terms in F2, see Eq. (BTT i. can result in a 
disentanglement over a finite time interval. 

Fig. |5]displays the concurrence using different initial con- 
ditions, by playing with the parameter a. This modifies the 
weight of the zero and double excitation states in the dynam- 
ics. When a = and K = the entanglement is maximal 
(C = 1) at all times. For finite K, keeping a = 0, it dies at 
the particular point at which |(ps)oi,io| = 0. Beyond that, it 
recovers to a value close to 1 . When we include the Q states, 
e.g., by taking a = 0.2, we observe the effect of entanglement 
sudden death, over a certain time interval. The duration of this 



interval grows when a is further increased up to a < 1/2. Be- 
yond this point the coherence in the P subspace may dominate 
over the population in the Q subspace, resulting in a positive 
value for F\ , eliminating entanglement sudden death. The be- 
havior at intermediate-strong system-bath coupling, K = 0.6, 
is included in Fig. |6j demonstrating that temporal oscillations 
are washed out. The dynamics at even larger K is similar in 
trends, with reduced concurrence value. 

The role of the temperature is displayed in Fig. [7] At high 
temperature the concurrence is zero. At intermediate values, 
T < 6 ~ 1 we find that its sole effect is a shift down of 
the qubit entanglement with increasing temperature. All other 
features (birth time, oscillation) stay intact. The simulation 
could not be performed at temperatures below T ~ 0.1 due to 
convergence issues in QUAPI. 

We can readily study the concurrence under different initial 
conditions for the P and Q subspaces, not necessarily in the 
form of Bell states, as long as Eq. ( TT~4b is obeyed. In particular, 
using a diagonal state for the time-zero reduced density ma- 
trix, similar features as those discussed above were obtained. 
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FIG. 5: Concurrence between the two qubits as a function of time, 
using Bell states [Eq. J22I with a — (full), a = 0.2 (dashed), a = 
0.5 (dashed-dotted) and a = 0.8 (dotted). Left Panel: K = 0.05. 
Right Panel: K = 0. Other parameters are the same as in Fig. [2] 
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FIG. 6: Same as Fig.|3]but at strong system-bath coupling K — 0.6, 
a = (full), a = 0.2 (dashed), a = 0.5 (dashed-dotted) and a = 0.8 
(dotted). 



V. UNIVERSAL FEATURES AT LONG TIME 

The long time behavior of the concurrence, representing the 
equilibrium limit, is displayed in Fig. [8] as a function of both 
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FIG. 7: The role of the bath temperature on the concurrence evolu- 
tion. T — 0.1 (dashed line), T = 0.2 (dashed-dotted line), T = 0.4 
(dotted line) and T = 0.6 (full line). We use Bell states [Eq. <22l l 
with a — 0.5 and K — 0.1. Other parameters are as in Fig. [2] 



K, the system-bath coupling parameter, and the initial state 
preparation ratio a, see Eq. d22l >. We note that the concur- 
rence can be significant in both the weak and strong coupling 
regimes, as long as the system evolves predominantly in either 
the P or Q subspaces. We now show that at weak coupling, 
K <C 1, and at low temperatures, T < Jj, for a broad range 
of parameters (as we explain below), the following general 
result holds 



C a<i(t 



1 - 2a. 



(26) 



The important implication of this result is that to the lowest 
order in K the concurrence deviates from unity due to the oc- 
cupation of the zero and doubly excited states in the system. 
This trend was observed before, e.g., in Ref. l29ll . However, 
here, for the first time, it is justified analytically, based on the 
spin-boson model behavior if3~3~ll . We derive Eq. d26*l > by study- 
ing the long-time limit of F\, as it dictates the concurrence 
when a < |, see Eq. (l24l . In the biased case, weak cou- 
pling theory (beyond the noninteracting blip approximation) 
provides Bill 



(Qz) = (As)oo,oo - (As)i 



a-r— tanh ( — - 



(27) 



A 2 



A 



<■/./■ 



is 



in the thermodynamic limit. Here A 2 = e 2 
a nontrivial function of K, cj c , and the bare tunneling element 
in the Q subspace, Jj. In the weak coupling limit we can 
write A e f / ~ J7, thus A& ~ y / e 2 + J 2 7 2 . 
polarization, we obtain the relevant term, 



Manipulating the 



\J (as)oo,oo(as)ii,ii ~ "J 1 - (j^J 



lanlr ).(28) 



The other element in F\ is the coherence in the P subspace. 
In the long time limit it satisfies 1I31I1 



l(As)o 



1,10 



i — a j /n 

tanh — 

2 Q \T 



(29) 



where il 2 — J 2 + 2J 2 K[i\ the proportionality factor obeys 
/x =_tb(iJ / nT) — hx(J/T) with ip as the digamma func- 
tion 113111 . As expected, the equilibrium concurrence depends 



on the environmental temperature, leading to entanglement 
degradation at high T, as observed in Fig. [7] Considering the 
low temperature case, T < J, J7, we note that the trigono- 
metric term in both Eq. (l28l l and Eq. d29l is close to unity. If 
we further work in the region e < Jj, the square root expres- 
sion in Eq. ((28) gets close to 1 . Under these broad conditions, 
the concurrence reduces to 



C.<4(* 



3) ~ (l-o) 



1 



1 - 2a - fiK(l - a) 



(30) 



One should note that the K dependence is more subtle than 
the simple linear scale attained here, since the tunneling ele- 
ment J should be corrected by K in a nontrivial manner 13111 . 
The simple result ( f30b provides us with some basic-interesting 
rules for building a long-time concurrence within the range of 
parameters mentioned above: (i) It decays linearly with the 
overall population placed in the Q (zero and double excita- 
tion) subspace. (ii) The reservoir temperature does not signif- 
icantly affect it. (iii) It does not depend on the qubit interac- 
tion energy. We note again that these observations are valid 
for a < 0.5, as long as T < J, Jj, e < J7 and K <C 1. 
When a > 0.5, the concurrence is determined by the compe- 
tition between F± and F2, see Eq. d20b . The numerics then 
suggests that C a> i(t — > 00) ~ 2a — 1 holds, for similar 
energy parameters. 

We conclude this section by emphasizing the implication 
of Eq. ( l26| i: One could set the steady-state entanglement in a 
dissipative system, by controlling the initial population in the 
P and Q subspaces. 




FIG. 8: The long time concurrence representing equilibrium behav- 
ior, for different initial states and system-bath coupling parameters. 
T = 0.2, J = 1, 7 = 0.5, 8 = 0.1 and ej = e 2 = 0.2. The long 
time limit was taken here as t — 100. 



VI. CONCLUSIONS 

Using exact numerical tools, we simulated the time evolu- 
tion of two qubits immersed in thermal environments, consid- 
ering a class of initial states for the subsystem. This task was 
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achieved by reducing the two qubits-bath model into two spin- 
bath systems, whose dynamics could be readily followed sep- 
arately. Using Wootters' formula for the concurrence lfl6ll . we 
quantified the degree of the qubits entanglement in time, ex- 
posing rich dynamics, including oscillations, delayed sudden 
birth, sudden death, and revival. Specifically, we showed that 
the occurrence of entanglement sudden death can be traced 
down to the initial population of the zero and double exci- 
tation states. The steady-state behavior was discussed in the 
weak coupling limit. 

Our results are significant for several reasons. First, we 
exposed a general mapping between an interacting two-qubit 
system embedded in a bath, and two spin-bath models, allow- 
ing us to simulate the dynamics of the original model using a 
numerically exact method that was developed for studying the 
prominent spin-boson case. Second, based on our mapping 
scheme, we calculated the concurrence measure and demon- 
strated the essential role of the environment in generating a 
stationary entanglement between the (interacting) qubits. By 
engineering the environment and its interaction with the sys- 
tem one could tune the degree of disentanglement in the sys- 
tem lioll . Earlier studies in this field have mostly treated a 



simpler version of our model, ignoring qubit-qubit interaction 
energy, further utilizing perturbative treatments. To the best 
of our knowledge, our work is the first to calculate the con- 
currence exactly in a dissipative and interacting qubit model. 

Future studies will focus on the dynamics of non-classical 
correlations beyond the entanglement measure, evaluating 
quantum discord ll32tl . This could be done by relying on the 
Y-form of the reduced density matrix H33l 13411 . With this at 
hand, we plan to study the dynamics of classical and quan- 
tum correlations in the qubit system, specifically, to investi- 
gate classical and quantum decoherence mechanisms and the 
possible transition between them j35ll . 
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